i) Missing data exploration (Figure not included in the article)

• Total number of missing data

# Missing data
#-------------

#Identify where data are missing
miss_dat <- trait_mat %>%
  # Remove the Bioturb_R modality that is not express by any of our polychaetes
  select(-Bioturb_R) %>%
  select(-Kingdom,-Phylum,-Class,-Order,-Genus,-Species) %>%
  bind_shadow() %>% 
  select(Family,species,ends_with("_NA")) %>%
  gather(variable,value,-Family,-species) %>%
  filter(!variable %in% c("Family_NA","species_NA")) %>% 
  mutate(variable=str_replace(variable,"_NA","")) %>%
  mutate(variable=str_replace(variable,"[.]","-")) %>%
  left_join(.,trait_list,by=c("variable"="Accronyme")) %>%
  mutate_if(is.character,as.factor) %>%
  mutate(variable=factor(variable, levels=unique(variable), ordered=T),Trait=factor(Trait, levels=unique(trait_list$Trait), ordered=T))

# Number of missing data
summary(miss_dat$value) %>%
  as.data.frame() %>% 
  set_colnames("Number of data") %>%
  rownames_to_column("Status") %>%
  mutate(Status= recode(Status, "NA" = "Missing", "!NA" = "Available")) %>%
  kable("html") %>%
  kable_styling(full_width = F, position="left")

• Modalities with missing data

trait_mat %>%
  select(-Kingdom,-Phylum,-Class,-Order,-Family,-Genus,-Species,-species) %>%
  miss_var_summary(.) %>%
  filter(pct_miss > 0) %>%
  dplyr::rename("Modality" = "variable", "Number of missing data" = "n_miss", "% of missing data" = "pct_miss") %>%
  kable("html") %>%
  kable_styling(full_width = F, position="left")

• 50% of missing data for Life span \(\Rightarrow\) Trait removed from the analysis
• 9% of missing data for Reproduction frequency \(\Rightarrow\) Missing data were imputed

• Distribution of missing data

p <- ggplot(miss_dat,aes(x=variable,y=species,fill=value))
p <- p + geom_tile(alpha=0.8)
p <- p + scale_fill_manual(values=c("#1a9641","#d7191c"), labels=c("Available","Missing"),name="Data availability")
p <- p + theme_light() + theme(text=element_text(size=24),axis.text.x=element_text(angle=90,hjust = 1),strip.text.x=element_text(angle=90),strip.text.y=element_text(angle=0),strip.text=element_text(size=26,face="bold"))
p <- p + facet_grid(Family~Trait,scales="free",space="free") + xlab("")

plot(p)

ii) Format the fuzzy-coded trait matrix

# Format the fuzzy-coded trait matrix
#-------------------------------------
trait_fuz <- trait_mat %>% 
  # Remove the empty bioturbation modality and the life span trait with lots of missing data
  select(-Bioturb_R,-ends_with("life_span")) %>%
  # Select only the species in the community matrix
  filter(species %in% colnames(com_pol)) %>%
  # Order the species as in the community matrix
  slice(match(colnames(com_pol),species)) %>%
  # Define row names and format the final trait matrix
  as.data.frame() %>%
  set_rownames(.$species) %>%
  select(-Kingdom,-Phylum,-Class,-Order,-Family,-Genus,-Species,-species)

# Recode the matrix to give the same weights to all traits 
#----------------------------------------------------------

# Create the function to convert the fuzzy coded matrix (going from 0 to 4 for each modality) to a relative coding in which all modality of a given trait sum to 1 for each species
transfo_traits <- function(x) {
  if (is.na(sum(x))) { # If we have NA's, gives back NA's
    rep(NA, length(x))
  } else if (sum(x) == 0) { #  If we have only 0's, gives back only 0's
    numeric(length(x))
  } else { # Else, we get the modality's score and divide it by the sum of the scores of the species for all modalities of this trait
    x / sum(x)
  }
}

# Prepare the trait list to apply the function
trait_list <- trait_list %>%
  # Remove the empty bioturbation modality and the life span trait with lots of missing data
  filter(Accronyme != "Bioturb_R" & Trait != "Life span") %>%
  # Reorder the list according to the order of the modalities in trait_fuz
  mutate(Accronyme=str_replace(Accronyme,"-",".")) %>%
  slice(match(colnames(trait_fuz),Accronyme)) %>%
  # Reorder the trait factor levels
  mutate(Trait = factor(Trait, levels=unique(Trait))) %>%
  droplevels()

# Compute the number of traits
n_trait <- nlevels(trait_list$Trait)

# Define the start and end point of each trait
n_mod <- trait_list %>%
  # For that we compute the number of modalities for each trait
  group_by(Trait) %>%
  summarise(n=n()) %>%
  pull(n) %>%
  # And we do the cumulative sum
  cumsum(.) %>%
  c(0,.)

# Apply the function
trait_prop <- trait_fuz 
trait_prop[,] <- NA

for(i in 1:n_trait){
  # Select all the modalities of the trait i
  tmp <- trait_fuz[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)]
  # Apply the function to this trait
  tmp <- t(apply(tmp, 1, transfo_traits))
  tmp <- as.data.frame(tmp)
  # Save
  trait_prop[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)] <- tmp
}

summary(trait_prop) %>%
  kable("html") %>%
  kable_styling(full_width = F, position="left") %>%
  column_spec(column = n_mod+1, border_right = TRUE, include_thead = TRUE)

# Create a vector defining the weight to be given to each modality so that all traits have equal weights in the analyses
#-------------------------------------------------------------------------------------------

mod_w <- trait_list %>%   
  # Compute the number of modalities in each trait
  group_by(Trait) %>%
  summarise(n_mod=n()) %>%
  ungroup() %>%
  # Add this to the trait list
  left_join(.,trait_list) %>%
  # Weight each modality so that the sum of the weights equals to 1 for each trait
  mutate(mod_w = 1/n_mod)

# Vérification que les poids somment bien à 1 pour chaque trait
mod_w %>%group_by(Trait) %>% summarise(w_tot = sum(mod_w))

## # A tibble: 10 x 2
##    Trait                                    w_tot
##    <fct>                                    <dbl>
##  1 Maximum size (mm)                            1
##  2 Feeding method                               1
##  3 Food size                                    1
##  4 "Preferred substrate position (adults) "     1
##  5 Living habit                                 1
##  6 Adult movement capacity (daily)              1
##  7 Bioturbation                                 1
##  8 "Sexual differenciation "                    1
##  9 Development mode                             1
## 10 Reproduction frequency                       1

# Format the weight vector
mod_w_tot <- mod_w %>%
  pull(mod_w) %>%
  set_names(mod_w$Accronyme)

iii) Missing data imputation

For the reproduction frequency, development mode, and sexual differentiation, data were missing for 9% (21 species), 7% (17), and 1% (3) of the species respectively. For these traits, we imputed missing values using nearest neighbour imputation relying on Gower dissimilarity that accommodates missing data. The imputation approach used in this article was originally developped for continuous traits (which is applicable here with our fuzzy-coded matrix) and was found to be particurlarly suitable when the trait distribution was unbalanced (Taugourdeau et al. 2014). Missing traits were imputed based on the median value of the functionally closest species for which the trait was known as well as those falling within a threshold dissimilarity of \(0.01\times d_{gower}\), with \(d_{gower}\) the dissimilarity between this closest species and the species to be inferred. This procedure gave similar results to imputation based on the 5 nearest neighbours using the kNN function. The species used to infer each missing data were then verified by experts of benthic taxonomy to ensure the ecological soundness of this imputation procedure.

Gower’s dissimilarity was computed on the relative scores (and not the raw fuzzy-coded scores) to account for the relative strengths of affinities within each traits (an affinity of 2 for a modality accompanied with an affinity of 2 for another modality of the same trait - meaning no preference between the two modalities of the trait - did not have the same meaning in our coding scheme that the same affinity of 2 accompanied with an affinity of 3 for the other modality - meaning a slight preference of the species for the second modality). Using the raw fuzzy-coded scores would not have accounted for this difference.

• Imputation procedure

# Compute the initial functional dissimilarity among species
dist_sp <- as.matrix(gowdis(x = trait_prop, w = mod_w_tot))

# Format the fuzzy-coded matrix
trait_fuz_init <- trait_fuz %>%
  rownames_to_column("species") %>%
  gather(Modality, value, -species) %>%
  left_join(.,trait_list,by=c("Modality"="Accronyme")) %>%
  select(-Modalities,-Type)

# Initiliaze the loop with the dataframe to be filled (_estim to remember that some trait data were estimated)
trait_fuz_estim <- trait_fuz_init

# We identify the species and traits for which data are missing (note that in this code, when data are missing for one modality, they are missing for all other modalities of the trait) 
trait_fuz_na <- trait_fuz_init %>%
  filter(is.na(value)) %>%
  distinct(species,Trait,value) %>%
  arrange(species)

# Identify the target species
target <- unique(trait_fuz_na$species)

# Object to record all modification that are done
logfile <- NULL

# Missing traits were imputed based on the median value of the functionally closest species for which the trait was known as well as those falling within a threshold of 1% of the dissimilarity between this closest species and the species to be inferred.
thresh <- 0.01

for ( i in target){
  
  # 1) For each target species, we identify which traits have missing data
  #------------------------------------------------------------------------

  trait_miss <- trait_fuz_na %>% filter(species==i) %>% distinct(Trait) %>% pull()
  
  # 2) We identify the functionnaly similar species
  #------------------------------------------------

  # Sort the distances of the target species with all other species
  sim_sp <- dist_sp[i,]
  
  # Remove the distance of the target species with itself
  sim_sp <- sim_sp[names(sim_sp)!=i]
  
  # Identify the closest species
  prox_sp <- sort(sim_sp)[1]
  
  # Identify the other species whose distance with the target species is within 1% of the distance of the closet species
  dist_sp_to_use <- sim_sp[ sim_sp <= prox_sp + prox_sp * thresh ]
  
  # Retrieve the names of the species whose value will be used
  sp_to_use <- names(dist_sp_to_use)
  
  # Do we have for these species at least one for which we have data for the missing traits 
  nodata <- trait_fuz_init %>%
    filter(species %in% sp_to_use & Trait %in% trait_miss) %>%
    group_by(Trait) %>%
    # If not, we put TRUE to absence, if we have, we put FALSE.
    summarise(Absence = all(is.na(value))) %>%
    pull(Absence)
  
  # If we have don't have any data for this trait for the functionnaly closest species and those around (within the threhsold), we print the table of the missing data for these species...
  if (any(nodata)){
    print("###########################################")
    print(paste("Missing data for the closest species that should have been used to estimate the missing trait data of species", i))
    print("------------------------------------------")
    print(trait_fuz_init %>% filter(species %in% sp_to_use & Trait %in% trait_miss) %>% select(-Trait) %>% spread(Modality, value))
  }
  
  # ... and we continue the procedure by looking at the following closest species (and those within 1% of its distance with the target species), and so on until we have data...
  while( any(nodata) ) {
    
    # Species already looked at
    n <- length(sp_to_use)
    
    # Identify the following closest species
    prox_sp <- sort(sim_sp)[1+n]
    
    # Identify the species within 1% of the distance of this "closest" species with the target species
    dist_sp_to_use <- sim_sp[ sim_sp <= prox_sp + prox_sp * thresh ]
    
    sp_to_use <- names(dist_sp_to_use)
    
    # Check data availability
    nodata <- trait_fuz_init %>%
      filter(species %in% sp_to_use & Trait %in% trait_miss) %>%
      group_by(Trait) %>%
      summarise(Absence = all(is.na(value))) %>%
      pull(Absence)
  }
  
  # 3) Get the median of their values (remove NA's if they are any)
  #--------------------------------------------------------------
  
  estim <- trait_fuz_init %>%
    filter(species %in% sp_to_use & Trait %in% trait_miss) %>%
    group_by(Trait, Modality) %>%
    summarise(estimate = median(value,na.rm=T))
      
  # 4) Replace the missing values (with the modalities to ensure everything is in the right order)
  #----------------------------------------------------------------------------------------
  trait_fuz_estim <- trait_fuz_estim %>%
    mutate(value=replace(value,species==i & Trait %in% estim$Trait, estim$estimate), Modality=replace(Modality,species==i & Trait %in% estim$Trait, estim$Modality))
    
 # 5) We save a record of the species used for the imputation, their distance to the target species and their trait values 
  #----------------------------------------------------------------------------------------  
  logfile <- trait_fuz_init %>%
    filter(species %in% sp_to_use & Trait %in% trait_miss) %>% 
    select(-Trait) %>%
    spread(Modality, value) %>%
    mutate(target=i, DISTANCE = dist_sp_to_use, Missing=paste(as.list(as.character(trait_miss)), collapse=" / ")) %>%
    select(target,species,DISTANCE,everything()) %>%
    bind_rows(logfile,.)
}

## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Asclerocheilus intermedius"
## [1] "------------------------------------------"
##                 species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Sclerocheilus minutus       NA         NA          NA          NA
##   Iteroparous Semelparous
## 1          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Goniadella spp."
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1  Labioleanira yhleni       NA         NA          NA          NA
## 2 Neoleanira tetragona       NA         NA          NA          NA
## 3    Pelogenia arenosa       NA         NA          NA          NA
##   Gonochoric Hermaphrodite Iteroparous Semelparous
## 1          4             0           4           0
## 2          4             0           4           0
## 3          4             0           4           0
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Neoleanira tetragona"
## [1] "------------------------------------------"
##             species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Pelogenia arenosa       NA         NA          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Nereiphylla rubiginosa"
## [1] "------------------------------------------"
##                species Iteroparous Semelparous
## 1 Notophyllum foliosum          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Notophyllum foliosum"
## [1] "------------------------------------------"
##                  species Iteroparous Semelparous
## 1 Nereiphylla rubiginosa          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Pelogenia arenosa"
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Neoleanira tetragona       NA         NA          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Pilargis verrucosa"
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Neoleanira tetragona       NA         NA          NA          NA
## 2    Pelogenia arenosa       NA         NA          NA          NA
##   Iteroparous Semelparous
## 1           4           0
## 2           4           0
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Polycirrus spp."
## [1] "------------------------------------------"
##                  species Iteroparous Semelparous
## 1 Amphitritides gracilis          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Scalibregma celticum"
## [1] "------------------------------------------"
##                      species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Asclerocheilus intermedius       NA         NA          NA          NA
## 2      Sclerocheilus minutus       NA         NA          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Sclerocheilus minutus"
## [1] "------------------------------------------"
##                      species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Asclerocheilus intermedius       NA         NA          NA          NA
##   Iteroparous Semelparous
## 1          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Sigalion mathildae"
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Neoleanira tetragona       NA         NA          NA          NA
## 2    Pelogenia arenosa       NA         NA          NA          NA

# Format the imputed matrix
trait_fuz_estim <- trait_fuz_estim %>%
  select(-Trait) %>%
  spread(Modality,value) %>%
  remove_rownames(.) %>%
  column_to_rownames("species") %>%
  # on réordonne les colonnes comme dans trait_fuz
  select(colnames(trait_fuz))

# Summary of the functional distances between the species to impute and those used to get the missing data
logfile <- logfile %>% 
  arrange(target, DISTANCE) %>%
  mutate_if(is.character, as.factor)

summary(logfile$DISTANCE)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.02679 0.03150 0.03699 0.04677 0.13439

The mean distance of the species used to estimate missing data is 0.04. At most, the imputation was made with a species with a dissimilarity of 0.13 with the species to impute. All these functional distances are fairly small compared to those observed across all the dataset (see below). In addition, they are in line with the \(d_{gower} \leq 0.05\) criteria used by Taugourdeau et al. (2014)

hist(dist_sp, main = "", xlab = "Gower distance between species pairs")
abline(v=summary(logfile$DISTANCE)[c(1,4,6)], col="red")
tmp <- summary(logfile$DISTANCE)[c(1,4,6)]
tmp[1] <- tmp[1] - 0.015
text(x=tmp, y=c(11000,11000,11000), labels=names(summary(logfile$DISTANCE)[c(1,4,6)]), col="red",xpd=T)

• Comparison with kNN

# Comparison with kNN imputation
#-------------------------------

# Imputation with median, 5 nearest-neighbours
trait_fuz_estim_kNN <- kNN(trait_fuz, weights = mod_w_tot, k = 5, numFun = median, weightDist = FALSE, imp_var = FALSE)

# Format the output
trait_fuz_estim_kNN <- trait_fuz_estim_kNN %>%
  mutate(species = rownames(trait_fuz)) %>%
  gather(variable,kNN_imputation, -species)

# Compare the two imputation approaches
compar_imputation <- trait_fuz_estim %>%
  rownames_to_column("species") %>%
  gather(variable,article_imputation,-species) %>%
  full_join(.,trait_fuz_estim_kNN)

# Differences between the two imputations
diff_imputation <- compar_imputation %>% filter(article_imputation != kNN_imputation)

# Number of missing data originnaly
n_miss <- trait_prop %>% gather(variable,value) %>% filter(is.na(value)) %>% nrow()

# Number of different imputations between the two approaches
paste("Number of different imputations between our approach and kNN:"," ",
nrow(diff_imputation),
"/",
n_miss,
sep="")

## [1] "Number of different imputations between our approach and kNN: 24/116"

# List of the species and modalities that would have change with kNN imputation
diff_imputation %>%
  mutate(diff = article_imputation - kNN_imputation) %>%
  arrange(diff) %>%
  select(Species = species, Modality = variable, Article_imputation = article_imputation, kNN_imputation) %>%
  kable("html") %>%
  kable_styling(full_width = F, position="left")

• Format the imputed matrix

# Recode the imputed fuzzy-coded matrix to give the same weights to all traits 
#----------------------------------------------------------

# Apply the function
trait_prop_estim <- trait_fuz_estim 
trait_prop_estim[,] <- NA

for(i in 1:n_trait){
  # Select all the modalities of the trait i
  tmp <- trait_fuz_estim[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)]
  # Apply the function to this trait
  tmp <- t(apply(tmp, 1, transfo_traits))
  tmp <- as.data.frame(tmp)
  # Save
  trait_prop_estim[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)] <- tmp
}

• Compare the trait-dissimilarity matrix computed from the initial and imputed trait matrix

dist_init <- gowdis(x = trait_prop, w = mod_w_tot)

dist_estim <- gowdis(x = trait_prop_estim, w = mod_w_tot)

qplot(dist_init,dist_estim) + 
  geom_pointdensity(alpha=0.5) + 
  geom_abline(intercept=0, slope = 1, col="red") +
  scale_color_viridis() +
  ylab("Dissimilarity from the \n imputed trait matrix") +
  xlab("Dissimilarity from the initial trait matrix")

iv) Assemblage‐by‐trait matrix

There are two ways to compute the assemblage-by-trait matrix containing the abundance of each modality in each site (also referred to as absolute mean value of modality, Fig.4 of Beauchard et al. 2017) :

1. One can first compute the \(CWM\) (Community Weighted Mean) matrix through the functcomp function of the FD package for instance. The CWM matrix represents the abundance-weighted mean affinity of a site to the trait modalities, computed as \(Mean~affinity = \frac{\sum Affinity_{[i]} \times Abondance_{[i]}}{Abondance_{Tot}}\). The latter can then be multiplied by the total abundance in each site to have the abundance of the modalities.

2. A more straightforward approach would be to do the matrix product of the site-by-species matrix with the species-by-trait matrix.

Both are be equivalent when there are no NA’s, which is the case given that we imputed all missing values. However, note that both approach would differ in the presence of NA’s. In the first approach, the CWM exclude NA’s in its computation. When multiplied back with the total abundance of each site (and if including in this total abundance the species for which trait data are missing), all traits will have the same total abundance (and missing data would then be considered as mean values). In contrast, with the second approach, the total abundance of the traits (and therefore their potential weights in some multivariate analyses) would vary according to their number of missing values (NA’s are just discarded). Note also that the functcomp function consider modalities with only 0 and 1 as binary by defaults. To treat these variables as numeric (and therefore take them into account in the CWM), we need to use the bin.num argument.

# Approach 1 : with the CWM
#-------------------------

# Compute CWM (we need to tell functcomp that trait_prop should considered as numerical) 
CWM_tot <- functcomp(x=trait_prop_estim,a=as.matrix(com_pol),bin.num=c("Bioturb_N")) %>%
  # Ensure that all variables are considered as numerical
  mutate_all(function(x){as.numeric(as.character(x))})

# Compute the abundance matrix
com_trait_estim_CWM <- rowSums(com_pol) * CWM_tot
rownames(com_trait_estim_CWM) <- rownames(com_pol)

# Check that the sum of abundances for each trait is equal to the total abundance of the sites
check_CWM <- com_trait_estim_CWM %>%
  rownames_to_column("id") %>%
  gather(Modality, value,-id) %>%
  left_join(., trait_list, by=c("Modality"="Accronyme")) %>%
  group_by(id,Trait) %>%
  summarise(value = sum(value)) %>%
  left_join(., data.frame(id=rownames(com_pol),tot=rowSums(com_pol))) %>%
  mutate(match_init = near(value,tot))

all(check_CWM$match_init)

## [1] TRUE

# Approach 2 : matrix product
#-----------------------------

com_trait_estim_MP <- as.matrix(com_pol) %*% as.matrix(trait_prop_estim)

# Equivalence between both?
#--------------------------

all.equal(as.matrix(com_trait_estim_CWM),com_trait_estim_MP)

## [1] TRUE

# Rename the matrix used in the following sections
#-------------------------------------------------

com_trait_estim <- com_trait_estim_MP

rm(com_trait_estim_MP,com_trait_estim_CWM)